Thermal Conductivity and Chiral Critical Point 
in Heavy Ion Collisions 

Joseph I. Kapusta^ and Juan M. Torres-Rincon^ 

^School of Physics & Astronomy, University of Minnesota, Minneapolis, MN 55455, USA 
^Departamento de Fisica Teorica I, Universidad Complutense de Madrid, 28040 Madrid, Spain 

(Dated: September 4, 2012) 

Background: Quantum Chromodynamics is expected to have a phase transition 
in the same static universahty class as the 3D Ising model and the liquid-gas phase 
transition. The properties of the equation of state, the transport coefficients, and 
especially the location of the critical point are under intense theoretical investiga- 
tion. Some experiments are underway, and many more are planned, at high energy 
heavy ion accelerators. Purpose: Develop a model of the thermal conductivity, 
which diverges at the critical point, and use it to study the impact of hydrodynamic 
fluctuations on observables in high energy heavy ion collisions. Methods: We apply 
mode coupling theory, together with a previously developed model of the free en- 
ergy that incorporates the critical exponents and amplitudes, to construct a model 
of the thermal conductivity in the vicinity of the critical point. The effect of the 
thermal conductivity on correlation functions in heavy ion collisions is studied in a 
boost invariant hydrodynamic model via fluctuations, or noise. Results: We find 
that the closer a thermodynamic trajectory comes to the critical point the greater is 
the magnitude of the fluctuations in thermodynamic variables and in the 2-particle 
correlation functions in momentum space. Conclusions: It may be possible to dis- 
cern the existence of a critical point, its location, and thermodynamic and transport 
properties near to it in heavy ion collisions using the methods developed here. 
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I. INTRODUCTION 

It has now been firmly established, via lattice calculations, that QCD with its physical 
quark masses does not exhibit a phase transition at finite temperature and zero baryon 
density but only a very rapid crossover from quark-gluon plasma-like behavior to hadronic 
gas-like behavior [l]-[5]. However, since the up and down quark masses are so small, and 
chiral symmetry is nearly exact, there are reasons to suspect that there is a curve of first order 
phase transition in the temperature T versus baryon chemical potential /i plane, terminating 
at a critical point at T^, > and /ic > 0. The existence of such a critical point has been 
found in various effective field theory models, such as the Nambu-Jona-Lasinio model [S]-[H], 
a composite operator model [9], a random matrix model [10], a linear a model |H], an effective 
potential model [H], and a hadronic bootstrap model [12]. Lattice QCD has yet to confirm 
or deny the existence of a critical point. The reason is that inclusion of a chemical potential 
does not allow for straight-forward Monte Carlo samplings of the field configurations. For 
reviews see Refs. and P^ . 

High energy heavy ion collisions may provide experimental evidence for a critical point 
and provide information on the behavior of the equation of state in its vicinity. Relevant to 
this are low energy runs at the Relativistic Heavy Ion Collider (RHIC), and in the future at 
the Facility for Antiproton and Ion Research (FAIR), at the SPS Heavy Ion and Neutrino 
Experiment (SHINE), and at Nuclotron-based Ion Collider Facility (NICA). Critical points 
are characterized by large fiuctuations. This led to the suggestion to study fiuctuations 
in conserved quantities, such as electric charge, baryon number, and strangeness on an 
event-by-event basis [I5l [16]. The effect is proportional to the spatial size of the domain 
or correlated region, which is probably rather small due to the finite size and lifetime in 
heavy ion collisions [17J. Therefore, it was suggested to measure higher moments to search 
for non-Gaussian behavior However, at present there is no experimental evidence for 
anomalous fiuctuations of this kind. 

A crucial issue is whether the critical point can ever be reached in a heavy ion collision. 
Colliding nuclei is necessary to reach high baryon density, but at the same time it creates 
entropy. If the initial entropy per baryon is much larger than that at the critical point, 
then the expanding matter will never pass close to it, even under the assumption of an ideal 
adiabatic expansion. The problem is similar to that of trying to create superheavy nuclei 
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by colliding nuclei: too much entropy is created. 

In this paper we construct a semi-realistic model for the thermal conductivity due to 
an assumed critical point in the QCD phase diagram. At a critical point the thermal 
conductivity diverges, as does the shear viscosity. However, the critical exponent for the 
shear viscosity is much smaller than that for the thermal conductivity, with the implication 
that the influence of the divergence of the shear viscosity is confined to a very narrow window 
in temperature, probably too small to have any effect on the matter produced in a heavy 
ion collision. In contrast, the temperature window for the enhancement of the thermal 
conductivity is much wider. These transport coefficients, along with the bulk viscosity, 
control the strength of hydrodynamic fluctuations in heavy ion collisions [H] . In particular, 
the strength of the hydrodynamic fluctuations due to thermal conductivity A are quantified 
by the correlation function 



where f{x) is a dimensionless fluctuation, or noise term, that appears in the hydrodynamic 
equations. Reference [19] applied the relativistic theory of hydrodynamic fluctuations to 
heavy ion collisions, with a specific example worked out for matter created with zero average 
baryon density.^ In this paper we will focus on thermal conductivity and ignore viscosities 
for simplicity of exposition. This allows us to study the influence of a critical point on the 
produced matter in a controlled and quantitative manner, although the model we use is not 
realistic enough for direct comparison with experiment. 

Reference [22] constructed the free energy in the vicinity of the critical point that in- 
corporated both the critical exponents and critical amplitudes. Using this free energy, the 
Landau theory of fluctuations was applied to estimate the probability of fluctuations away 
from the equilibrium state and were found to be very large. However, the Landau theory 
considers baryon number fluctuations in a finite volume in contact with a heat reservoir at 
fixed temperature, which does not adequately represent the space-time evolution of matter 
in a heavy ion collision. Dynamical simulations of spinodal decomposition, or phase separa- 
tion, were done in p3| and |24J, but without the incorporation of intrinsic hydrodynamical 
fluctuations or noise. In this article we assume that the entropy created in the collision is 

^ Earlier studies for extracting the shear viscosity in heavy ion coUisions were done by |20) and for the bulk 
viscosity by |21) . 
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too high to allow any trajectories in the T — jj, plane to pass through the coexistence region. 
Rather, following the earlier remarks on entropy in heavy ion collisions, we assume that the 
trajectories always pass to the left of the critical point, never entering the coexistence region 
or crossing the line of first order phase transition. This is a more conservative scenario for 
heavy ion experiments. 

Hydro dynamic fluctuations, or noise, may be crucial for studying the effects of a critical 
point in heavy ion collisions. As an analogy, Ref. [25j performed a theoretical study of 
the breakup of liquid nanojets with the conclusion that "noise is the driving force behind 
pinching, speeding up the breakup to make surface tension irrelevant" . Similar conclusions 
were reached in Ref. [26] , which studied the breakup of liquid nanobridges. 

The outline of the article is as follows. In Sec. 2 we construct a semi-realistic model for 
the critical enhancement of the thermal conductivity. This is based on the mode-coupling 
theory which has been successfully applied to the liquid-gas phase transition in a variety of 
ordinary atomic systems. In Sec. 3 we discuss the equation of state to be used in conjunction 
with the thermal conductivity to model the expansion phase of heavy ion collisions. In Sec. 
4 we implement these ideas in a boost-invariant hydrodynamic model; although not the 
most relevant fluid dynamic model for studying the critical point, it allows us to gain insight 
and intuition, and to discover the magnitude of the effect on observables in Sec. 5. We 
summarize and conclude in Sec. 6. 

It should be acknowledged that there are other sources of fluctuations in heavy ion 
collisions, such as initial state fluctuations, fluctuations induced by jets and other high 
momentum-transfer processes, and fluctuations during hadronization in the flnal state. 
These were surveyed in Ref. [19]. Fluctuations caused by passage near a critical point 
should be characterized by a strong beam energy dependence. Fluctuations due to jets 
should not be relevant at the energies where the critical point might be reached. 

II. THERMAL CONDUCTIVITY 

In ordinary materials, thermal conductivity is typically inferred from measurements of 
a diffusion constant rather than measured directly. There are several kinds of diffusion 
constants depending on the experimental conditions. The baryon diffusion constant -D^ is 
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defined by the diffusion equation 

f = O.V^X . (2) 

where X is the departure from the equihbrium baryon density n or baryon chemical potential 
II at fixed temperature T. The fully relativistic expression for this diffusion constant is 

AT /n\2 



(dn/dii)T 

where w is the enthalpy density and A is the thermal conductivity (dimension of energy- 
squared). The partial derivative in this expression is just the baryon number susceptibility 

XB = idn/df,)T . (4) 

The isothermal compressibility kt is 

where P is the pressure. These are simply related by xb — ti^kt- 
The diffusion equation for heat is 

(6) 

which is carried out at constant pressure. This diffusion constant can be derived from first 
order relativistic viscous fiuid dynamics to be 

Dt^-, (7) 
Cp 

where cp = T{ds/dT)p is the heat capacity per unit volume at constant pressure. 

For liquids and gases near their critical point it is easier to measure Dt than Db- As 
is well-known, the thermal conductivity diverges at such a critical point. It is conventional 
and useful to separate the thermal conductivity into a smooth background piece A^ and a 
piece AA which diverges at the critical point and which goes to zero away from the critical 
point so that A = A* + AA. Exactly how this is done is not unique and is considered a bit 
of an art. Then AA = cpADx- 

We would like to know the thermal conductivity not only in the asymptotic critical region 
but in the non-critical region too. For this one must go beyond the renormalization group, 
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which is only vahd asymptotically close to the critical point. A rather successful approach 
to this problem is to use mode-coupling theory for the dynamics of critical fluctuations in 
fluids extended away from the critical region The basic idea is the recognition 

that in a fluid the slow modes are the diffusive modes of heat and viscosity while the sound 
mode is considered a fast mode. The approach is quite general, and has been developed 
for fluids near and above the critical temperature in the universality class of the liquid-gas 
phase transition. Only the essential points will be reviewed and summarized here [32] . 
The part AZ^t due to the critical enhancement takes the form 

= ^^^-^^ • 

Here ^ is the correlation length, t] is the shear viscosity, Rd is a universal constant ap- 
proximately equal to 1.05, and is a crossover function which goes to zero as qoC, goes to 
zero and which goes to 1 as qoC, goes to infinity. The qd is a cutoff in wavenumber and is 
material-dependent^. As the critical point is approached, ^ — )■ oo, and the Stokes-Einstein 
diffusion law is recovered 

ADr^f^. (9) 

The smooth background value of rj, without the critical enhancement, is sometimes used 
because its critical exponent of 0.063 for the reduced temperature [33]- [25] is much smaller 
than the critical exponent for the thermal conductivity, and its divergence at the critical 
point contributes negligibly to the results. 

A model for the crossover function Qlq^C,), which partly accounts for non-asymptotic crit- 
ical behavior for the thermal conductivity, was presented in Fig. 4 of |32j- The complicated 
numerical results can be approximated to two-significant digits in the range 0.5 < qo^ < oo 
by 

1 04 

n{x) ^ 0.48tanh(0.23x) + arctan(0.65x) . (10) 

TT 

The correlation length needs a precise definition. Ref. [32] uses 

(U) 



^ To be precise, both AD^ and the crossover function fl depend on the wavenumber q. In the hydrodynamic 
hmit one takes — > while qo^ remains finite. Hence we consider ^{qo^) a function with a single 
argument. 
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where 



Xb 



Pr 



Xb 



dn\ 



(12) 



is a dimensionless susceptibility. Both and F are smoothly and slowly varying functions 
of n and T. A function 



T, 



ref 



rcf J 



T 



(13) 



was used in ^32j so that the enhancement from criticality goes to zero at some reference 
temperature Tref. The critical exponents are 7 ~ 1.24 for the isothermal compressibility 
or baryon number susceptibility, and u ~ 0.63 for the correlation length. These are the 
critical exponents for systems in the same static universality class as the liquid-gas phase 
transition and the 3D Ising model; it is generally accepted that the QCD phase transition 
under discussion belongs to the same universality class [3S]- It has further been shown [3Zj 
that this QCD phase transition belongs to the dynamic universality class of model H of Ref. 
[3H]. 

An explicit expression is needed for the correlation length. Based on the work in [22] we 
take the baryon number susceptibility in the critical region to be 



^X*i 



5Pr 



This results in the correlation length 



6-1 
2^ 



An 



+ 55 \r] 



<5-l 



-1 



(14) 



e(n,T)=6 



6 



7 



An 



f + 55 \r] 



5-1 



-I//7 



(15) 



Here t = {T — Tc) /Tc > and 7] = {n — nc)/nc (not to be confused with the shear viscosity). 
The other critical exponent is 5 = 4.815. The An is the discontinuity in the baryon density 
at zero temperature. It was estimated in [22j as An = ric/S. All the other constants or 
slowly varying quantities are absorbed into 



eo = Co 



5Pr 



(5 + i)/.r 



(16) 



Note that the expression (14) automatically goes to zero, or becomes very small, far away 



from the critical point. Hence there is no absolute need to make a subtraction as in (13). 
Finally, when considering heavy ion collisions later on, it is necessary to know AA for t 
somewhat less than zero but still outside the phase coexistence region. To do so we shall 
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replace V in the previous formulas with which is a smooth extrapolation outside the 
domain of the derivation. 

Now we come to a discussion of the remaining parameters. Ref. |32] made a theoretical 
fit to the transport properties of carbon dioxide and ethane. Their fit parameters are: 
^0 = l-SA, = 2.0A, and T = 0.0481 for carbon dioxide, and ^0 = 1-9A, q^^ = 1.83A, 
and r = 0.0541 for ethane. (Note that light scattering can be used to independently fix the 
critical amplitude ^0 in these materials, but such is not possible with the matter produced in 
high energy heavy ion collisions.) The first point to notice is that ^0 ~ Qd^ for both tances. 
They are smaller than the average particle separation of ric^^^ = 5.39A(carbon dioxide) 
and ric^^^ = 6.23A(ethane), and comparable to the physical dimensions of the molecules. 
We will take qo^o = 1 and assume that the cutoff is qo = ttTq where Tq is the crossover 
temperature at zero chemical potential; see the next section. For the other parameters we 
choose r = 0.05 and from [22] = bPc- Finally, with Tq = 170 MeV (for example) we get 
^0 = 0.37 fm and ^0 = 0.69 fm. 

There are various theoretical approaches for calculating the background thermal con- 
ductivity A''. A calculation in QCD to lowest order in as but to all orders in In as gives 
[39] 

0.165 1 

~ «2in(o.497/«,)T ■ ^ ^ 

This is for two quark flavors; similar results were obtained for three flavors. This result 
requires that ^ 1 which is valid only at asymptotically high temperatures. No calcula- 
tions have been done with nonzero chemical potentials with this accuracy. Earlier, Ref. [10] 
estimated the thermal conductivity with a chemical potential for two flavors as 

A'«.^(-)^ (18) 



a^ln(l/as) \n, 

using the relaxation time approximation to the Boltzmann equation. 

The relaxation time approximation is especially useful for calculating the background 
thermal conductivity in the hadronic phase. From [?T] and [IQ] 

A' ^ ^X:(2.. + l)e'.^/-/ (ff - . (19) 

where the sum is over all particles i with baryon number 6j, spin Sj, and single particle 
energy Ei = ^/p'^ + mf. The Tj(p) is the momentum-dependent relaxation time. In the limit 
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that the baryon chemical potential goes to zero this results in 




Hence the baryon number fluctuations, as reflected by ([T]), remain even when the average 
baryon density is zero. Note that in this limit only baryons contribute to the thermal 
conductivity. This limit will be a good approximation when the enthalpy per baryon is 
much greater than the nucleon mass. 

It should be remarked that in some special situations the number of particles, in particular 
pions, is held fixed due to relatively slow particle creation and annihilation reactions. In 
these situations the net baryon density is replaced by the total particle density, for example 
ra^. See, for example, [H] and [12] • In the context of unitarized chiral perturbation theory 
the thermal and electrical conductivities, which are related through the Wiedemann-Franz 
law, have been calculated in |13]. In [H] the electrical conductivity has been related to the 
total photon yield udN^/d^p in the limit pt — )■ 0. However, due to our higher temperatures 
-where pion number is not conserved- we will not follow that path here and consider net 
baryon density. 



III. EQUATION OF STATE 

We need a smooth background equation of state for two reasons. First, we need the 
isobaric specific heat in order to convert the thermal diffusion constant to the thermal con- 
ductivity in the critical region. Second, we need an equation of state to solve the background 
hydrodynamic equations to find the space-time evolution of the independent thermodynamic 
variables T and /i. 

The isobaric specific heat diverges at the critical point. It obeys the thermodynamic 
relation 

The isochoric specific heat has critical exponent a = 0.11 while the susceptibility has critical 
exponent 7 = 1.24. Hence the critical exponent for the isobaric specific heat is 7 and for the 
thermal conductivity it is 7 — = 0.61. This is under the reasonable assumption that the 
coefficient of xb in the expression for cp — cy has a nonzero finite value at the critical point. 
It should be noted that use of the equation of state which includes the critical amplitudes 



(21) 
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and exponents as given in [22] results in the vanishing oi cp — cy- In particular, the function 
/o(t) in that paper is only given to first order in T but it is needed to second order to 
compute cp. This just means that a smooth background equation of state must be used to 



calculate the coefficient of xb in (21). 

In |22] it was assumed that the critical point lies somewhere along a crossover curve 
parameterized by 

with an estimate that Tq ~ 170 MeV and /io ~ 1200 MeV. The high energy density equation 
of state was taken to of the form 

P = A^T^ + AsTV^ + Ao^^^ - CT"^ - B . (23) 

This equation of state was used to calculate the critical density, critical entropy density, and 
so on once was given. 

If either P or e is assumed to be constant along the crossover curve then one has the 
following condition on fio/To. 

^0^-^4 + ^4 = 0. (24) 

There is also a condition on C. If P is constant along the crossover, as argued in |22] then 

C = f^l (^2 - 2Ao^) . (25) 

For definiteness we shall take the coefficients Ai for those of a noninteracting gas of massless 
gluons and Nf flavors of massless quarks. 

TT^ / 21iVf~ 

18 ' 

Ao = . (26) 

Using Nf = 2 and To = 170 MeV results in /io = 1218.48 MeV. At the critical point 

/ T^2\ -1/2 





and 



^=(^)'^ (28) 
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FIG. 1: The phase diagram showing the crossover curve and the zero pressure curve. 



for the both this equation of state and for the parameterization near the critical point 
constructed in [;22]. Finally, as in |22] we take B = O.STq. The phase diagram is illustrated 
in Fig. [T] 



Using (23) the terms appearing in (21) are 



9/i 
df 



-2A2Tfi 



s 

n 



T /2A4T2 + A2/i2-C 



(29) 



(30) 



These are smooth in the vicinity of the critical point and provide a nonvanishing coefficient 
of xb- Note that if there was no term in the pressure, C = 0, then an adiabatic expansion 
would imply that T//i = constant. In that case the factor 

n 2 



9/i 
df 



+ 



n 



12 





FIG. 2: Plots of AA using the parameters given in the text. 

is constant during adiabatic expansion and cooling. 

Combining the results from the previous section and this one, we have found a represen- 
tation for AA. In Fig. [2] we plot it as a function of 1] for fixed values of t and in a contour 
plot as a function of t] and t. Note that the energy scale is GeV, not tens or hundreds of 
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MeV. Therefore one should expect significant fluctuations if the trajectory of the expanding 
matter takes it near the critical point. In the remaining part of this article we will neglect 
the background thermal conductivity and use only AA. If the magnitudes of the resulting 
correlations are observably large, then they would only be greater if was included. 

To conclude this section we should point out that one can add higher derivative terms to 
the dissipative part of the baryon current in the Landau-Lifshitz frame of reference. This is 
discussed in Appendix [X} 

IV. FLUCTUATIONS IN BOOST INVARIANT HYDRODYNAMICS 

In this section we will study the effects of fluctuations on the expanding and cooling 
matter produced in heavy ion collisions if the trajectory in the T — fi plane passes near the 
critical point. For this purpose we will use the 1+1 dimensional boost-invariant (Bjorken) 
hydrodynamic model, similar to what was done in Ref. [19]. This model is too idealized to 
make any sort of direct comparison to experimental data. In addition, this model assumes 
highly relativistic beam energies, probably much higher than is required to achieve the 
high baryon densities necessary to probe the critical point. Nevertheless, it does provide 
some guidance and intuition before one attempts to study the problem with much more 
sophisticated and numerically intensive 3+1 dimensional viscous fluid dynamics. Since the 
method and equations are so similar to those in [T9]| we will leave out some of the details. 
The gaps can easily be filled in with a little effort. 

The energy-momentum tensor in ideal fluid dynamics is 

T'"' = wu^u" - Pg'"' . (31) 

The shear and bulk viscosities are ignored to focus on the effects of thermal conductivity. 
In the boost-invariant hydrodynamics the flow velocity has the nonvanishing components 

M° = cosh(^ + oj) , 

= sinh(^ + u) . (32) 

Here ^ is the space-time rapidity (not to be confused with the correlation length) and u{^, r) 
is a fluctuation that depends on both ^ and the proper time r. The baryon current is 

J/^ = nu^" + /'^ , (33) 
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FIG. 3: Evolution of baryon density (left) and entropy density (right) as a function of time. They 
fall as l/r so that the entropy per baryon is constant within the hydrodynamic expansion. 

where is a fluctuation, as described in [19]. The smooth, background fluid equations lead 
to the simple equations of motion 

rl Q e 

(34) 



ds s ^ 
— + - = 
ar T 



and 



dn n 
— + -=0 
dr T 



(35) 



independent of the specific equation of state. The solutions are 



s(r) = SiTi/r 



(36) 



and 



n(r) = UiTi/T 



(37) 



where Sj and rii are the entropy and baryon densities at some initial time Xj. Some represen- 
tative solutions for s and n for later use are represented in Fig. [3j The initial temperature 
is Tj = 250 MeV, the initial time is Tj = 0.5 fm/c, and the initial chemical potential was 
chosen to be fii = 420, 620, and 820 MeV for the three cases. The adiabatic trajectories, 
corresponding to these three cases, and are also shown in the phase diagram in Fig. |4} 
Trajectories I and II represent crossover transitions while trajectory III passes very close to 
the critical end point, which is here chosen to be located at = 160 MeV and fic = 411.74 
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FIG. 4: The phase diagram showing the crossover curve and the three trajectories used in the 
computation. 

McV. The entropy per baryon at the critical point is 19.96, while for trajectories I, II and 
III it is 37.98, 26.08 and 20.06, respectively. The time evolution is stopped when the zero 
pressure curve is reached, which is at t/ = 3.04, 3.30, and 3.68 fm/c, respectively. In reality, 
matching to a full hadronic equation of state should be done, but we don't do it for this 
illustrative example. 

The full equations are linearized in the fluctuations, such as 5n, and these in turn are 
linear functionals of The nonvanishing components are 



I' = s(T)/(e,T)sinhe, 
I' - s(r)/(e,r)coshe 



(38) 



on account of the condition that u^I^ — 0. Notice that / is dimensionless as the entropy 
density has been factorized out for convenience. The average value of /(^, r) is zero and the 
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correlation with itself was given in Eq. Q. The linearized equations are 

, 



dSe _ duj 
T— + dw + w— 
or 



dSn did df 

T— + 6n + n— + s— = , 



(39) 
(40) 

Here n and w are the smooth background solutions which depend only on r. Note that the 
noise drives the baryon fluctuations, and if the average baryon density is zero there is no 



dr 

d{wuj) 



d5P 

+ 2WUJ + ^ir:- = . 



coupling to u and so there is only one equation (40) to solve. 



It is useful to identify the independent variables in Eqs. (39)-(41 ) as 5s, 5n, and co. Then 
they take the form 



dSs 
dSn 



. du fisdf 



dtu df 
or a4 a4 

r— + 1 - vl)uj + — ^ + — ^ = . 



(42) 
(43) 
(44) 



Here is the physical, adiabatic speed of sound squared, while = {dP/de)n and 
{dP/de)s. They are related by 

2 ^ Tsvl + finvl 
w 



(45) 



See Appendix |B} Finally, it is convenient to use dimensionless independent variables. The 
linearized equations of motion are written in terms of 6s /s, 6n/s and u: 

d f 5s\ dtu fi df 



^ dr \ s 

d f 5n 
dr \ s 



+ 



d^ Td^ 
n du df 



s dt, dt, 



du 



vlTs d fds 



w d^ 



+ 



Ins d ( bn\ 



- =0 



w d^\s) 



In Fourier space 



/oo 
rfee-^'^«X(e,r) 
-oo 



for any variable X. Then equations (46)-(48) can be expressed as a Langevin equation: 

+ Dt^ + /ii = , 



(46) 

(47) 
(48) 

(49) 
(50) 
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where the vector ip, the drift matrix D and the stochastic noise term n read: 

I — ^ / n n ik \ 



s 

6h 



s 

\u J 



D 



ikv, 















^ — 

n „ , 


ikVg — 



n 



ik- 



l-vl 



n 



ik 



T 
1 





(51) 



V ^ / 



WW ' 

The homogeneous equation is solved by the evolution operator U(A;; r, r') that is calculated 

as 



\}{k\ T, t') = Texp 



" dr" 



Bik,T"] 



(52) 



where T is the time-ordering operator. Once that is known, the solution to the inhomoge- 
neous equation can be expressed as 



ip{k,T) 



dr' ~ 

—V{k;T,T')n{k,T')f{k,T') 

T 



(53) 



with the assumption that the solution at the initial time is zero. Up to this point no 
specific equation of state has been assumed and the results are completely general (within 
the context of the hydrodynamic model). 

In general, the solution to the above equations cannot be written down in closed form. 
As a consequence, the evolution operator generally would need to be computed numerically 
once an equation of state is specified. The exception is when the drift matrix is constant in 
time. That was the situation studied in [19], and it would be the situation here too if the 



equation of state is given by (23) with C = 0. If the drift matrix was time-independent, 



T){k,T") = D(fc), then the solution would be 

lJ{k-r,T') = 



D(fc) 



(54) 



We will approximate U by the above expression with D evaluated at r. This should be 
a very good approximation at large T and/or /i, and at least semi-quantitatively valid in 
any case. (The reason for evaluating the drift matrix at r instead of r' is discussed more 
later.) To avoid this approximation would require taking C = 0, and it is arguable whether 
it would be a better approximation to the real physics. 

The calculation is further simplified as the drift matrix is diagonalizable. The eigenvalues 
are: 



do 
d± 



, 

a±/3 



(55) 
(56) 
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with 



The evolution operator reads: 



lJ{k- T, t') = \Jo{k- r') + ^p^lJ_{k; t') 



where 



ci+ — d. 
( 



jr'/rf 
d^ — d. 



-U+(^-;r') , 



(57) 



(58) 



Uo(A;;t') 



n/jLVg —pisVg ^ 



-Tnv^ Tsvl 



v 











(59) 



U±(fc;TO 



V 



v"^ Ts 

w 


-2—d^ 
w 


—ik 


vl Tn 

W 


w 


—ik— 
s 


—ikvl — 

w 


—ikvl — 
w 


—d± 



(60) 



The solution to the inhomogcneous Langevin equation gives the response functions 
Gx{k,T,T'). They are calculated from 



^ G,{k;T,T') ^ 

Gn{k;T, t') 
yG^{k;T, t') J 



U(A;;t, T')n{k,T') 



(61) 



to be 

Gs{k;T,T') 

Gn{k;T,T') 



ik n 



ik 



+ ( 



r 



sinh (j3 In + cosh In 



— sinh (j3 In — j + cosh (j3 In — j 



G.{k;r,r') = ~{vl-vl) 



,\ a 



sinh ( /3 In — j . 



(62) 



T / \ T' 

where, unless otherwise indicated, the variables on the right hand side are functions of r. 

Although 5n/s and 5s /s were the natural dimensionless variables to use when solving the 
equations of motion, it is very useful to look at the variables 



5P 



and 



^0- 



s J 



/ Sn\ n f 5s 
\s) s\s 



(63) 



(64) 
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The corresponding response functions are 



Gp{k]T,T') =ik^ {v'^s-vl) 



T 



a 



T 



T 



sinh ( /3 In — j + cosh ( /3 In — 



r 



and 



Gcr = ik 



w 



(65) 



(66) 



The reason that these are interesting and relevant is that a change in pressure at constant 
entropy per baryon corresponds to a sound wave, while a change in entropy per baryon at 
constant pressure corresponds to diffusive heat flow, and so these are physically orthogonal 
variables. This should be apparent mathematically as well. Notice that Gp contains only 
cosh and sinh terms; in the large k limit these become sinusoidal oscillations. On the other 
hand G^ has no such terms; in coordinate space it is the derivative of a Dirac delta function. 



Approximating the matrix D as a constant, as in Eq. (54), can be examined in light 
of its diagonalized form. The approximation will be good if |91nf^/(91nr| ^ 1. This is 
a good approximation because it is only the presence of a term in the pressure that 
causes v^. to deviate from its asymptotic value of 1/3. For a static, uniform system no such 
approximation is necessary. See Appendix [Cj 

In the presence of the fluctuating forces, the solution to the equations of motion for 
X = 6s/s,6n/s,oj,6P/Ts,6a reads: 



X(fc,r) 



dr' ~ 

—Gxik,r,r')fik,T') 



(67) 



The correlation function of the fluctuating force ([T| is 



(/(ri,6)/(r2,6)) 
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5{xi - X2) 



An 
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and where A is the transverse area in the Bjorken model. The Fourier transform is 
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(68) 



(69) 
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In Fourier space the correlation function of the fluctuating variables then reads 



{X{h,n)Y{k2,T2)) 
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(71) 
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The equal-time correlator at the final or freeze out time r/ is 



CxY{k]Tf) 



~A 



dr 



A(r) 



n{r)T{T) 
s{t)w{t) 



GxY{k;Tf,T) 



where 



Gxvik; Tf, t) = Gx{k] rj, r)G'y(-A;; r/, r) 
Finally, the correlation function in space-time rapidity is 



A 



''f dr 



n{T)T{i 



SiTjWiT) 



(72) 



(73) 



(74) 



where GxY^ii — ^2;t/,t) is the inverse Fourier transform of GxY{k;Tf,T). As mentioned 
earlier, we will ignore the background thermal conductivity A'' in our numerical results. The 
final correlations due to thermal conduction would only be greater if A^ was included. 

As discussed in some detail in [19], the response functions generally contain step functions, 
Dirac delta functions, and first and second derivatives of Dirac delta functions. These are 
due to the fact that signals created by the fiuctuating sources propagate at the speed of 
sound; the singular parts are located at the sound horizon in space-time rapidity. The 
distance a signal can propagate in space-time rapidity between time ri and T2 is 

dr 



r 



(75) 



Two points separated by twice this distance would receive a signal at the same time if it was 
sent from a source midway between them. As in one can subtract the singular terms 
to study the regular behavior of this function. 

GxUk\ r/, r) = Gxy(fc; r^, r) - G'^^{k- r^, r) . (76) 

We show some examples of G'^y(^; r/, Xj) in Fig. [s] These represent the wake behind 
the propagating sound front. The singular expansions of these coefficients are given in 



Appendix D Note that the subtraction of the singular part is only made for illustration 
purposes. The final particle correlation contains the entire GxY{k',Tf,T) function. It is 
noteworthy that if we approximate the drift matrix at the time r' instead of at r there 



would be a small singular piece in Gpp at the distance given by (75) in addition to at twice 



that distance. This is unphysical, and reaffirms the choice of the time r rather than r'. 
It is interesting to observe that there are no sound waves emitted if Va = Vg = Vn, only 



pure diffusion, as may be seen from Eq. (62 ). Such is the case for the simple equation of state 
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P = v'^e — B where is a constant. Sound waves are generated only if the pressure responds 
differently to variations in energy density depending on whether the entropy density, the 
baryon density, or the entropy per baryon is held fixed. It is immediately apparent that 
QCD does not have the property that all three speeds are equal. The equation of state used 



for the space-time evolution in this article, given by Eq. (23), has different speeds because 
it contains a term in the pressure. For very large energy densities it has a decreasing 
influence, and all three speeds approach 1 / a/3; see Appendix tel 



V. PHENOMENOLOGY 

In this section we consider observable consequences of fluctuations caused by the existence 
of a critical point in the QCD phase diagram. Although direct comparison to experiment is 
not to be expected, we will find that the effects are quite significant and definitely worthy 
of further study. 

The analysis follows that described in [19] very closely. For this reason, we only give the 
key steps here. The number of particles with degeneracy d per phase-space volume is 

^=(2^^^"'P^' 

with 

/(x,p)=e-(-^-'^)/^ (78) 
being the Boltzmann distribution function. The four-velocity of the fluid cell is 

= (cosh(^ + u), u_L, sinh(^ + u)) . (79) 

The distribution of particles at the freeze-out surface S j is given by the Cooper- Frye formula 

In the hydrodynamic model being used here 

d^cT/i = Tf d^ d'^x±m± cosh(?/ - ^) . (81) 
The variable y represents the particle rapidity 

p'^ = (m_L cosh y, p_|_, m_L sinh?/) , (82) 
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where 

m± = \Jrn? + p\_ (83) 

is the transverse mass. The number of particles per unit rapidity is then 

dN _ dAvf 



4) 



^^^y J d^cosh{y-^) j rf^p_Lm_Lexp{-[(cosh(?/-^-w)m_L -/i)/T]} , (8 

where the integration over gives the transverse area of the coUision A. If we neglect all 
fluctuations {5T = 6(1 = u = 0) we get the average of dN/dy as 

{^) ^ J^y2^^^^ J cosh(?/-^) J dp_Lp±m±exp{-[cosh{y - ^)m±]/T} . (85) 
In order to perform the integration over p_i_ we use the following formulas: 

r I 1 

/ dp^p^m^e-^""^ = - e-'^"'[2 + 2cm + [cmf] = —1(3, cm) , (86) 
J 

[ dp^p^mle-'""^ = \ e-"'"[6 + 6cm + 3(cm)2 + {cmf] = ^1(4, cm) . (87) 
J 

At the freeze-out time we obtain 

dN_ ^ cL4^ „,/T, r ^f. m 



. o e^^^ ^ / 3, — cosh a; . (88) 

dy 47r2 J _^ cosh^ x \ 'Tf ) ^ ' 

Now we consider fluctuations of dN/dy and eventually its two-point correlation. To do so, 



we expand the exponential term in (84) to first order in fiuctuations around the freeze out 
value: 

T = Tj + 5TiTj,0 , (89) 
/X = fif + 6fi{Tf,0 . (90) 

The Boltzmann factor becomes 

exp {-[{cosh{y - ^ - u)m± - /i)/T]} exp {-[(cosh(y - ^)m± - Ai/)/T/]} 

X |l + ^[m^cosh(y-0-^/]+u;(e)^sinh(y-0 + ^| , (91) 

where the fiuctuations are understood to be evaluated at Tf. The fiuctuation in the number 
of particles per unit rapidity is then 



^ /'dN\ dA 



dy J (27r)3 



Tf J d^cosh{y-^) j d^p^m^ exp {-[(cosh(y - ^)m_L - /i/)/T^]} 
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X 



f mo 



I -^2- ["^^ cosh{y - - ^/] + w(0 ^ sinh(?/ - + J • 



(92) 



We need to replace the fluctuations 6T and 5/i by 6s and 5n. We use 

ST = ^Ss-^Sn, 
A A 



XT^l r- . Xtt ^ 
OS -\ — on 



(93) 



A A 

where xtt = d^P{T, ij,)/dT^, xt^ = d'^P{T, ^j)/dTd^ji, and x^/. = d'^P{T, ij,)/diJ,'^, and where 



(94) 



We now perform the integration over p_L with the help of (86) and (87). The fluctuation of 
dN/dy reads: 

'dN\ dATfT] 



' dy ) Atx'^ 



^ e^-f/^-f / d^ 



-Fs{y - + ^F^{y - + —Fn{y - 

s s 



Here we have introduced the functions 

Fs{x) ^ -^20^T (A,^coshx] - ^^^^+^y^/^^ r f 3,^cosha: 
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A cosh^ X 
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Finally, we construct the rapidity correlator: 



dATfTj 
47r2 



2\ 2 



e^'^^/^^- 1 d^, I rf6 



X 



X,Y 



where the sum runs over 6s/ s, 6n/s, and u, and where 



(95) 

(96) 
(97) 

(98) 



(99) 



Cxy(6-6;r/) = (X(a;r;)r(6;r;)) 



(100) 



To simplify calculation of the 2-particle correlation function it is convenient to use the 
Fourier transforms. 



dyij \dy2 



dATfT: 



'2\ 2 



47r2 



111] ^^f^f/Tf 



I ^,^''^yJ2Fx{k)Fy{-k)CxY{k;Tj) 
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(101) 
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with Ay = Hi — y2- The expression for Cxy is: 

Cxvik) = -j -A(r) J Gx{k; Tf, T)Gyi-k; tj, t) . (102) 

To ehminate the transverse area we normahze the correlation function by the event average 
of dN/dy to obtain 

2 



where 



C{Ay)= fdke''^yY,Fx{k)Fy{-k) r^x(—) Gx{k;Tf,r)Gyi-k;Tf,T) (104) 
and 

N {m/Tf) ^ r -^T U ^ cosh x) . (105) 
J —oQ cosn X \ J J 

Now we are prepared to study the influence of the hydrodynamic fluctuations on the 
2-particle correlation functions. Figure [6] shows the thermal conductivity as a function of 
time for the three trajectories chosen earlier. As the trajectory passes closer to the critical 
point the enhanced thermal conductivity is probed ever more closely. The effect of this on 
the 2-particle correlation functions is shown in Fig. [Tjfor protons [d = 2, nip = 939 MeV). 
For charged pions (d = 2, m^r = 138 MeV, zero chemical potential) we show in Fig. |8]the 2- 
particle correlation function with no chemical potential fluctuation {6fi = 0) (left panels) and 
with pion chemical potential fluctuation equal to that of the protons (right panels). Which 
one is closer to reality cannot be determined without explicitly introducing a conserved 
electric current on the same footing as the baryon current. The presence of additional 
fluctuations enhances the magnitude of the correlation between particles. Both pions and 
protons exhibit the influence of the thermal conductivity. The shape is approximately the 
same for all trajectories. The protons have a maximum at Ay = and a minimum near 
Ay = 0.95; the pions have a maximum at Ay = and a minimum near Ay = 1.7 (with 
chemical potential fluctuations) or near 1.3 (without chemical potential fluctuations). The 
magnitude of the correlation, on the other hand, increases dramatically as one goes from 
trajectory I to II to III. As mentioned earlier, the boost invariant hydrodynamical model 
is not sufficiently realistic to consider any comparison to experiment, especially since one 
needs a large initial baryon number. Nevertheless, the size of the effect in this model bodes 
well for future theoretical studies and experiments. 
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FIG. 8: Particle correlation function for pions at zero chemical potential. Top: Non-zero fluctua- 
tions of chemical potential. Bottom: No fluctuations of chemical potential. 

VI. SUMMARY AND CONCLUSIONS 

We applied mode coupling theory, together with a parameterization of the equation of 
state that incorporates the correct critical exponents and amplitudes, to develop a model 
for the thermal conductivity in the vicinity of the critical point. This contribution to the 
thermal conductivity incorporates the correct critical behavior but can be used in the non- 
asymptotic region as well. The thermal conductivity quantifies the strength of particular 



28 



hydro dynamic fluctuations via the fluctuation-dissipation theorem. To obtain insight into 
what effects might result in heavy ion collisions as a consequence of the critical point, we 
studied a simple boost-invariant hydrodynamic model. We conservatively assumed that the 
entropy per baryon created in a heavy ion collision is too large for an adiabat to enter 
the mixed phase. The thermal conductivity along the adiabatic trajectory is enhanced the 
closer the trajectory comes to the critical point. The flyby of the critical point results in 
fluctuations in the temperature, baryon chemical potential, and local flow velocity which 
evolve with time and are not the same as fluctuations in the initial conditions. We found 
that the growth of the thermal conductivity near the critical point implies the existence 
of two-particle correlations over 2 units of rapidity for protons. Moreover, the strength 
of this correlation increases the closer the expansion trajectory comes to the critical point. 
With the inclusion of other transport coefficients (shear and bulk viscosities) this correlation 
can be further enhanced. We have found that the fluctuations in the baryon density - 
in particular the chemical potential fluctuations- are the most important effect for the 
magnitude of correlations compared to the thermal and velocity fluctuations alone. For this 
reason, we think that the usually neglected baryon diffusion coefficient -alternatively, the 
thermal conductivity- could be of interest for accessing the critical behavior by two-particle 
correlations. However, their practical use relies on the ability of the heavy-ion factories to 
produce trajectories that pass close to the critical point. 

There are many natural extensions to this study. For example, inclusion of the regular 
part of the thermal conductivity would increase the magnitude of the fluctuations and hence 
correlation functions in the flnal state observables. To quantify fluctuations due to electric 
charge, pions for example, requires the introduction of the electric charge current in addi- 
tion to the baryon current. Introduction of the strangeness degree of freedom -the kaon 
multiplicity is typically larger than proton multiplicity- would also enhance the correlation 
function as one includes additional fluctuations in the strange chemical potential. Inclusion 
of the thermal conductivity into the fluid equations of motion, not only in the fluctuations, 
should certainly be done. This could potentially increase correlations because it couples 
some of the fluctuations among them in such a way that the drift matrix does not contain 
any nonvanishing element. The effect of this modiflcation for the shear and bulk viscosities 
has been studied in [IS] leading to a smoothing of the singularities appearing at the sound 
horizon. 
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An extension to the region of first order phase transition would be interesting. Decreasing 
the initial entropy per baryon below its value at the critical point would result in instabilities 
and phase separation, with probably greater consequences than a flyby of the critical point 
as studied here. The importance of the thermal conductivity in this regard has already been 
noted |15]. Ultimately a 3+1 dimensional calculation for the fireball evolution is necessary 
to study the hydrodynamical correlations not only in rapidity but also as a function of the 
azimuthal angle, analogous to what was done in Ref. |16]-|18] for initial state fluctuations. 
Their effects on the flow harmonics and the angular correlation of particles ought to be 
important. Such a study is of course much more difficult than the one presented here since 
one must perform extensive numerical calculations. 
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Appendix A: Baryon Current in a Gradient Expansion 

In the Landau-Lifshitz definition of flow velocity the baryon current is modified by a term 
linear in a gradient and proportional to the thermal conductivity. 

.r = nu^' + \ A^(/3^) , (Al) 

where /3 = 1/T and A'^ = — u^{u ■ d). The tensorial structure is determined by the 
requirement that the dissipative term be orthogonal to the flow velocity so that m • J = 0. 
Going to second order in a gradient expansion yields 

= nu^ + \(^^ [(1 - tbm ■ 9)/3/i] , (A2) 

where tb is a new time constant which is in principle dependent on T and /i. For baryon 
diffusion at fixed temperature and in the local rest frame this results in a modified diffusion 
equation. 

— = DbV dfi - tbDbV . (A3) 

This is a causal equation, unlike the conventional diffusion equation. It should be mentioned 
that in the Bjorken hydrodynamic model that dissipative term in the baryon current van- 
ishes identically because A'^ is acting on a function of r only. The dissipative term would 
contribute to the linearized fluctuation equations, but the resulting fluctuations would be 
one higher order in A. 



Appendix B: Speed of Sound 



The adiabatic speed of sound (constant a = s/n) can be computed for (23). It is w ^ = 
I - Af ^ with 

2 2 CT a{A2T' + 6Aofi')-2A,fiT 

^"- = 9d(T:^ ^tt~^ ' ^^^^ 

where 

rf(T, /x) = 2A4A2T4 + (I2A4A0 - ADfi^T^ + 2^2^0/1^ - lCiA2T^ + QAofi^) . (B2) 

Notice that Av^ is proportional to C and hence vanishes for a conformally invariant equation 
of state. Since ds = adn, the speed at constant baryon density, v'^ = {dP/de)n, can be 
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obtained by taking cr — > oo. The speed at constant entropy density, Vg = {dP/de)s, can be 
obtained by taking cr — > 0. It is easily verified that 

w 

a relationship that is independent of the specific equation of state. Of course waves do 
not physically propagate at constant n or s, only at constant cr, but these definitions are 
useful for various intermediate steps in various applications. For example, a thermodynamic 
identity is 

dt.\ _vlcy-S ^ ^^^^ 



dT J ^ n 
so that 

cp - Cy = V^CyTxB ■ (B5) 



Appendix C: Perturbations of a Static Uniform System 

It is instructive to be reminded of perturbations of a static, uniform system at rest. The 
analog of the boost invariant hydrodynamics used as an example in the text is a three 
dimensional system with perturbations dependent on z and t but independent of x and 
y. The resulting equations for the fluctuations are solved straightforwardly since the drift 
matrix is time-independent. The response functions correspond directly to those given in 
Section IV with the following replacements: a — > 0, cosh — > cos, sinh sin, ln(r/r') — )■ t—t', 
and (3 — )■ v^^k. Also the function a;((^,r) is replaced by the velocity perturbation v{z,t) in 
the 2;-direction. 

Gs = -^i^{vl+{vl-vl)cos[kvAt-t')]} , 

Gn = '^{vl+{vl-vl)cOs[kv^{t-t')]} , 

~ Jv S 

Gv = {vl- vl) sm[kv^(t - t')] , 

Gp = ik^ {vl - vl) cos[kva{t - t')] , 

Ga = ik^^ . (CI) 

From these it is immediately clear that a pressure disturbance will travel with the speed of 
sound Va while disturbances in the entropy per baryon are local in coordinate space. 
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Appendix D: Singular Part of the Response Functions 

The response functions GxY{k',T,T') contain terms that produce singularities when the 
inverse Fourier transform is performed. The general expansion of the singular part is ob- 
tained by a Laurent expansion in and retaining the regular terms at /c = 0. For the 
cases XY = ss, nn, sn, uu, the general expansion of the singular part reads: 

G'^y'{k;T,T') = (aik^ + hi) + (a2fc' + ^2) cos{2v^Lk) + ^ sin(2^;^Lfc) 

K 

+ {a4,k^ + 64) cos{vaLk) H ^\ii{vaLk) , (Dl) 

with L = log(r/r'). The coefficients aj,6j are polynomials in L: 

fli = ^ aijU , &i = ^ hjU , (D2) 

i i 

whose coefficients Ojj and hij depend on some thermodynamical quantities. The expansions 
of Qi and hi are represented in Tables |T] and [IT] for the cases XY = sn and XY = uu. Note 
that for the case XY = ss and XY = nn one should use Table [T] and make D ^ A, F ^ 
C,E ^ B and A^D,C^F,B^E, respectively. The cases XY = suj,nuj follow a 





j = 


j = l 


i = 2 


i = 3 


aij 


-AD - CF/2 











hij 


-BE/2vl 











a2j 


-CF/2 











b2j 


BE/2vl 


{CE + BF)a^/2vl 


CFa^/4:vl 





azj 


-{CE + BF)/2v„ 


-CFo?l2v^ 








hj 


-{CE + BF)a'^/Avl 


-BEo?l2vl - CFa^/Svl 


{CE + BF)a'^/Avl 


CFa^/l2vl 


a^j 


-CD - AF 











hij 





{BD + AE)a'^/2vl 


{CD + AF)a^/Svl 





a^j 


-{BD + AE)/v„ 


-{CD + AF)a'^/2va 








hj 


-{BD + AE)a^/2v^ 


-{CD + AF)a^/^vl 


{BD + AE)a^/^vl 


{CD + AF)a^/4&vl 



TABLE I: Coefficients of the singular expansion (Dl) for the case XY = sn. 
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j = 


j = 1 


J = 2 


i = 3 


aij 


GV2vl 











hi 













a2j 


-GV2vl 











b2j 


-G^ay2vi 





G'^a^/Avt 





asj 





-G^ay2vl 








hj 





-bG'^a^/Svl 








TABLE II: Same as Table U for the case XY = ujuj. 



different singular expansion. They generically read: 



Gx^{k;T,T') = cik + C2k cos{2vaLk) + (csk'^ + sm(2vaLk) 
+ c^k cos{va^Lk) + (csA;^ + ^5) sm{vo-Lk) . 

The functions Cj, rfj admit an expansion in powers of L: 



(D3) 



(D4) 



In Table III we show the coefficients Cj,-, dij for the case XY = stu. Finally, the case XY = nuj 



uses the same expansion but with the change B^E,C^F, A^D. 





j = 


i = 1 


J = 2 


J = 3 


Clj 


iBG/2vl 











C2j 


-iBG/2vl 


-iCGo?l2vl 








C3j 


iCG/2v„ 











dsj 


iGGo?l^vl 


-iBGo?l2vl 


-iGGa^l^vl 










-iAGo?l2vl 










iAGjva 











d5j 


iAGo?l2vl 





-iAGa'l^vl 






TABLE III: Coefficients of the singular expansion (D3) for the case XY = su. 
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For the ease of simplicity, in this Appendix we have defined the following quantities: 

--^(^)". 0^li^-<)i^'- 



